/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		basic_geometry.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-03-23 15:33
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-07-20 12:46
*
*   Description:
*
=============================================================================*/
#ifndef BASIC_GEOMETRY_H
#define BASIC_GEOMETRY_H

#include "../types.hpp"
#include <numeric>

namespace geoxfem {
namespace geometry {

//(-pi,pi]
inline _real anti_clock_angle_from_x_axis(const _real& v1x,
                                          const _real& v1y,
                                          const _real& v2x,
                                          const _real& v2y) {
  return atan2(v2y - v1y, v2x - v1x);
}

//(-pi,pi]
template <typename point_type>
inline _real anti_clock_angle_from_x_axis(const point_type& v1,
                                          const point_type& v2) {
  return atan2(v2.y() - v1.y(), v2.x() - v1.x());
}

template <typename point_type>
int if_at_left_wing(const point_type& p0,
                    const point_type& p1,
                    const point_type& p) {
  _real x0 = p0.x();
  _real y0 = p0.y();
  _real x1 = p1.x();
  _real y1 = p1.y();
  _real x  = p.x();
  _real y  = p.y();

  _real side = (x1 - x0) * (y - y0) - (y1 - y0) * (x - x0);

  return (side > 0) ? 1 : -1;
}

template <typename point_type>
point_type getCentroid(const std::vector<point_type>& geometry) {
  std::vector<point_type> plist = geometry, _plist;
  std::sort(plist.begin(), plist.end());
  std::unique_copy(plist.begin(), plist.end(), back_inserter(_plist));
  return std::accumulate(_plist.begin(), _plist.end(), point_type{})
         / static_cast<_real>(_plist.size());
}

struct minus_pi_to_pi {};
struct zero_to_2pi {};

// calculate the 2D angle between p0-pc-p1
template <typename point_type>
struct anti_clock_angle {
  anti_clock_angle() = default;
  anti_clock_angle(const point_type& p0,
                   const point_type& p1,
                   const point_type& pc)
      : m_corner{pc},
        m_start_direction{anti_clock_angle_from_x_axis(pc, p0)},
        m_angle_radian{anti_clock_angle_from_x_axis(pc, p1)
                       - m_start_direction} {}

  template <typename RangePolicy>
  _real rad() const {
    throw "not implemented yet";
  }
  template <>
  _real rad<zero_to_2pi>() const {
    return m_angle_radian < 0 ? m_angle_radian + M_PI * 2 : m_angle_radian;
  }
  template <>
  _real rad<minus_pi_to_pi>() const {
    if(m_angle_radian > M_PI)
      return m_angle_radian - M_PI;
    if(m_angle_radian < -M_PI)
      return m_angle_radian + M_PI;
    return m_angle_radian;
  }

  template <typename RangePolicy>
  _real deg() const {
    return rad<RangePolicy>() * rad_to_deg;
  }

  bool ifObtuse() const {
    return abs(m_angle_radian) > M_PI_2;
  }
  bool ifAcute() const {
    return abs(m_angle_radian) < M_PI_2;
  }

private:
  point_type m_corner;
  _real      m_start_direction;
  _real      m_angle_radian;
};

}  // namespace geometry
}  // namespace geoxfem

#endif /* BASIC_GEOMETRY_H */
